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INTRODUCTION 



The method most frequently used for solving numerically the random resistor network 
(RRN) problem has changed over time surprisingly often: relaxation methods for solving 
Kirchhoff's equations were adopted in the seventies, while the early 80's were the time of 
the random walk method; then the transfer-matrix (TM) approach [|T| came into fashion; 
and next the node-elimination method came forth in the 90's 0,0]- A "Fourier acceleration" 
method was also proposed in mid-80's |Q. Renewed interest in direct methods to solve the 
set of Kirchhoff equations arose after the paper of Edwards et al. had been published in 
1988 1^: standard algebraic multigrid(AMG) method generally used for solving large linear 
sparse systems was applied. In a recent paper ^ the standard Kirchhoff's set was reduced 
by Green's function formulation of Kirchhof 's laws. 

The random walk method is probably the worst among the methods listed above. Al- 
though in some applications the random walk method could be more suitable than the others 
the main reason for its frequent use appears to be the nice exposition given in Stauffer's 
famous introductory book [0]. This method faces the same problem as many iterative meth- 
ods for solving the Kirchhoff's equations: its performance decreases rapidly at the critical 
region of a metal-insulator-like phase transition - s.c. critical slowing down(CSD). Random 
walkers diffuse anomalously slow at criticality(p = Pc), hence the diffusion constant (i.e. the 
conductivity) estimations require more computer time at p = Pc than for p > pc- In the 
same way network size scaling at criticality leads to a faster increasing of numerical efforts 
than the number of resistors involved. The origin of CSD is not so transparent when the 
Kirchhoff's equations are iteratively solved. CSD amounts here in increasing the number of 
iteration needed to reach a certain precision. Probably CSD stems from the fractal geometry 
of the resistor network. Such geometry leads to multifractal distribution of voltage drops 
across the net (if an external voltage is applied) and in this way reduces the speed of 
convergence in iterative solvings. 

The AMG, transfer-matrix and node elimination methods are free of CSD in a sense 
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specific for eacli metliod. Tlie speculation that AMG metliod eliminates the critical slow- 
ing down completely (or almost completely) relies mostly on numerical data [^] which is 
far from exhaustive. The node elimination method calculates directly (without any voltage 
evaluations) the network conductance by applying consecutively the star-triangle transfor- 
mation in a way to reduce the number of sites in the network until one resistor is left. The 
computational effort is proportional to the number of resistors so it is faster a.t p = pc than 
for any p > pc- In a similar manner the TM approach is faster at the critical point for a 
given system size, but here computations scale with size in a non-trivial way which will be 
discussed further. It is important to note that a modification of the TM approach in order to 
evaluate the voltage drops distributions cite is possible 0; such modification is impossible 
within the node elimination method. 

For all of these three methods, the "preconditioning" of the system by extracting the 
connected (spanning cluster) or bi-connected (percolation backbone) component could sig- 
nificantly improve the method performance. For node elimination and AMG approaches 
such extraction could be easily implemented. 

In this paper I am concentrating on the TM approach, presenting a modified algorithm, 
which allows preconditioning by extraction of spanning clusters. 

A conceptually important feature of any TM approach is that one does not have to 
consider the entire system (or its states) at a time in order to calculate its physical properties. 
Typically, one only requires information about state n in order to proceed to state n+1 and 
subsequently discards the information about state n. In contrast, the known ways [p|,p!0|-[T2 



to extract the backbone requires that the percolation structure is kept in computer memory 
in its entirety. 

I present a TM algorithm which is improved in comparison to the previous TM formula- 
tions in two ways. It has inherited the important feature of "volt age- source book-keeping" 
from an earlier modification of the "canonical" TM approach made by the author ||T3| for 
application to quasicrystalline and random lattices. This feature makes possible a better 
utilization of the dilute structure of the random networks. Second, the system is precondi- 



tioned by spanning clusters extraction. The specific way of extraction reduces significantly 
the memory requirements, which otherwise are very restrictive. 



I. THE TRANSFER-MATRIX APPROACH 

The TM approach to the numerical solution of the RRN problem has been presented 
first by Derrida & Vannimenus in 1982 ^ and has been elaborated subsequently by several 
groups [P^SH TTD. Characteristic for the TM approach in two dimensions (2D) is the use of 



infinitely long strips of finite width L cut from the resistor network and, analogously, in 
three dimensions (3D) the use of "bars." The similarity to the transfer-matrix method of 
the statistical mechanics of spin systems consists of the introduction of a matrix A{M) 
which represents the properties (in the RRN problem the conductivities — see below) of the 
semi-infinite strip between — oo and strip slice M. As, e.g., in a strip slice of a resistor 
network on the square lattice may consist of the vertical resistors in column M and the 
horizontal resistors which connect columns M and M — 1. Knowledge of the conductivity 
matrix A{M) and the resistor configuration in slice M + 1 is sufficient to calculate A{M + 1) 
for the next slice. Iteration for all subsequent slices finally obtains the conductivity of the 
whole strip. In the case of a resistor network of unit resistors and insulators the long edges 
of the strip are thought of as electrodes and the resistors in the upper and lower layers are 
taken to have zero resistance. 

If the resistor network has a fractal structure its conductance tends to zero as the system 
size increases — in analogy to its mass density decreasing to zero. More quantitatively, 
the infinite spanning cluster at the percolation threshold Pc is a fractal for which finite-size 
scaling theory shows that its conductivity should scale with the system size L as L~^^^ , 
where v is the percolation correlation-length exponent and t is the percolation transport 
exponent. 

The TM approach has been used first for obtaining precise estimates of t for percolation 
on the square and cubic lattices |T^,[T^. In Refs. the matrix A{M) is updated after 
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the addition of every single resistor (the program is pubhshed in instead of using matrix 
equations as in [Q. Thus, the calculations are simplified and accelerated. 

In order to define the matrix A{M), we attach voltage sources to the open ends of the 
resistors at the right end column of the growing semi-infinite strip. 

The matrix A{M) is defined by attaching to the sites of the current right end column 
of the semi-infinite strip voltage sources Vi, where i labels the row position in the strip and 
thus assumes values from 1 to L, the width of the strip. Since the voltage-current relations 
in the network are linear the current from any selected sources, say source j, is a linear 
function of the voltages Vi, 

i^=j:MM)v,. (1) 

i 

The relation defines the matrix elements Aij{M) of the TM A{M). From now on I will 
suppress the argument M when no confusion can arise. 

When a horizontal resistor R is added to row k then the matrix A changes to A' with 
matrix elements 



For infinite i?, a case that we encounter in insolator-resistor mixtures, Eq. (Q) simplifies to 

= K - (3) 



When we add a vertical resistor between two adjacent sites k and / of a new column then 
four matrix elements change, 

A'„i = Aki - 1/R, 
A\, = A, - l/R, 

(4) 

A'kk = + l/R, 
= An + II R. 

From Eq. it is clear that, in the limit M — > oo, e.g., the difference A{M)ll — A{Q)ll 
tends to the transverse conductance of the strip of width L. From analysis of the conductivity 



scaling of strips with different L one obtains the conductivity scahng exponent. For the 
percolation cluster at Pc this exponent equals the ratio t/u. 

The advantages of the TM approach have been described in the pioneering works and 
p!6| . Here, I would like to point out to the reader its main drawback: the size of the matrix 
A and the computational effort grow very fast with the strip width L. 

In particular, the size of the matrix grows as and for 2D and 3D respectively. If we 
consider a site percolation model in 2D then addition of every new column leads to adding 
an average of p^iL — 1) horizontal and p'^{L — 2) + 2p vertical resistors. Taking the size of 
the matrix A into account, we find that Eq. (|^) is applied oc times whereas only oc L 
operations are required for Eqs. (^. Thus, it is clear that for widths larger than 10 — 15 
lattice spacings more than 90% of the time is spent on calculating the relations (Q). In 3D 
the situation is even worse: the upper bound for the computational efforts scales as L^. But 
in fact this bound is overestimated: Ref. [|1^] points out that the computational effort scales 
as due to the fact that the matrix A is sparse, i.e., most of its elements equal zero. 

In the next section, I will describe a modification of the TM approach that overcomes 
these problems in parts. 

II. THE MODIFIED ALGORITHM 

The site-percolation case will be considered without loss of generality. The "voltage- 
sources book-keeping" procedure is described in the next subsection. Second(Sec.II.B.) 
comes the method for extracting the spanning clusters and at last(Sec.III.C.) I present the 
main steps in the complete algorithm. 

A. Conductivity calculations 

Let us reconsider Eq. (|1]) in the case of a general resistor network. Let some network 
nodes of an arbitrary resistor network be connected with external voltage sources V^). Then 
Ij [Eq. (|l|)] is the current from source j and the absolute values of the off-diagonal elements 
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of matrix Aij,i ^ j; represent the conductances between voltage sources i and j. The 
diagonal elements give the conductance between the respective source and the "ground" — 
the other sources set to zero voltage. In fact, we can likewise interpret the matrix elements 
Aij introduced in Eq. ([1|) in the previous section. The difference is that the number of 
sources in the present case does not strictly depend on the strip width and their connection 
to the right side end is not mandatory (see Fig. 1 in Ref. Il3|). Eqs. (0-^ still apply in the 



present general case if we allow in the general geometry the "addition of a vertical resistor" 
as connection with a resistor of two sites with voltage sources attached to them and the 
"addition of a horizontal resistor" as insertion of a resistor between a site and the voltage 
source previously attached to this site. In other words, "adding a horizontal resistor" creates 
a new site and moves the voltage source to it. 

Based on these general concepts, we may formulate as algorithm the conductivity calcu- 
lations Eqs. (0-^ for the resistor strip case and the construction of the strip. The algorithm 
has three main steps. 

(i) Add a new site to the right end of the already existing strip and attach a 
voltage source to this site. 

(ii) Find the neighbors of this site among the sites already present and connect 
them with resistors. Update the matrix using Eq. (Q). The algorithm should 
ensure that the neighbors have their own voltage sources attached, 

and (iii) if a site with its attached voltage source is located within the bulk of the 
growing strip then detach the source to free it for subsequent attachment to a 
new site on the growing edge of the strip. To this end, we (a) insert an insulating 
"resistor" between the network site and the previously attached source. Then 
(b) we update the TM according to Eq. (|^). The information about the voltage 
source index is kept on a stack, ready to be (c) used again when a new site is 
added to the right end of the strip. We have done nothing more than adding a 
bridge of infinite resistance between two points of the network which does not 



alter its conductivity properties and moved a voltage source. 

Since a regular lattice structure of the resistor bonds is not a prerequisite for the con- 
ductivity calculation outlined above, an algorithm based on the steps (i-iii) is useful for 
calculation of percolative conductivities of quasicrystalline and random lattices ||T3| , p^ . 



Since we have to only update the TM for lattice sites that are actually connected to the 
strip by resistors of finite resistance and since we always apply the simpler Eq. (j^) instead 
of Eq. (0), the outlined procedure is already faster than the standard algorithm p6[] . 

The matrix size, i.e., the number of matrix elements, instead of being is only ^ p^L^ 
for the square lattice and ~ p^L^ instead of L'^ for the cubic lattice. The scaling with L of 
the matrix size is not altered but the way of handling the voltage source numbers allow for a 
significantly smaller prefactor and (more important) it facilitates the system preconditioning 
which does lead to improvement of memory and performance scaling. 

B. Spanning clusters extraction 

We achieve an improvement of the scaling of the matrix size as a function of L by the 
extraction of spanning clusters. The spanning clusters are defined as the percolating clusters 
which connect the top and bottom edges of the strip or, respectively, the bottom and top 
faces of the bar in 3D. At the spanning clusters in strip geometries represent the incipient 
infinite percolation cluster. Its fractal dimension df is 91/48 in 2D and around 2.5 in 3D 
0. If only voltage sources connected to the spanning clusters contribute to the matrix size 
then this size should scale as L?'^'^i'~^\ where the exponent df — 1 reflects the system width 
dependence of the scaling of the spanning cluster's sites found in a — 1 dimensional cut 
p!9| . Thus, in 2D the matrix size scales as L^-'^^- instead of . In 3D the number of the 



matrix elements is oc L^ " rather than L^. 

How does one extract the spanning clusters in strip geometries? Several general algo- 
rithms exist to solve this problem [|T^-[T^. However, they all require that the percolation 



structure has been created beforehand and is stored in its entirety leading to large computer 
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memory requirements ( e.g. for a bar 10^ x 100 x 100 one has to consider 10^ sites and 10^ 
matrix elements). 

I now present an algorithm which partly resolves this problem which would otherwise 
limit the strip length. 

The method for extracting the spanning cluster is based on the Hoshen-Kopelman algo- 
rithm for cluster counting. As is well-known this algorithm requires only consecutive 



d — 1 dimensional cuts of the lattice to be kept during its lattice "sweeping." Cluster in- 
formation is stored in one ID array — the array of cluster sizes and pointers, sometimes 
denoted as the array of "labels of labels (LOL)" 0. The index into this array represents 
the cluster labels and its elements are either "cluster roots" — then containing the size 
of a specific cluster — or pointers to these cluster roots — i.e., negative numbers whose 
absolute value is equal to the index of the array element corresponding to the cluster root. 
Moreover, the cluster root element may be used to store other information about its cluster, 
e.g., whether the cluster touches the upper or/and lower layer of the strip. 

Running the Hoshen-Kopelman algorithm requires that the percolation structure be 
scanned twice in order to extract the spanning clusters. These two runs are required because 
if we reach a site during the first run we cannot decide if this site's cluster will eventually 
turn out to span. To avoid storage of the entire cluster in memory, we perform the second 
sweep based on a repetition of the pseudo-random number sequence that created the first 
sweep cluster. During the second sweep the LOL array is examined to decide which cluster 
a site belongs to and whether this cluster spans. Only sites belonging to spanning clusters 
enter the conductivity calculations. 

Thus, instead of storing the percolation cluster structure itself, we only store the LOL 
array. The key question for the proposed algorithm is the size of the LOL array that has to 
be retained in memory between the two Hoshen-Kopelman sweeps. To keep its size small 
I apply a procedure to recycle unused labels 1^1],^. The size of the resulting LOL arrays 
turns out to be less than 0.5% of the memory required to retain the whole cluster structure 
in memory. 
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c. 

The full algorithm, including both the spanning clusters extraction and the conductivity- 
calculation procedures, may be summarized as follows: 

1. scan the random structure with the Hoshen-Kopelman procedure constructing the 
LOL array by a label recycling technique. After the sweep, keep the LOL array in 
memory. 

2. repeat the scan using the same pseudo-random number sequence: 

after creation of a new site DO: 

decide by comparing the new and the stored LOL array whether 
this site belongs to a spanning cluster. If it does then the site 
enters the TM conductivity calculations. These proceed according 
to step (i-iii) as outlined in the previous subsection. 

3. When the second scan terminates calculate the transverse conductance per unit length 
as, 

^ A{M)ll - A{Mo) LL 



M-Mo 

where I have used Mq = M/5 to reduce boundary effects. 



(5) 



III. PERFORMANCE SCALING RESULTS 

I have developed the algorithm described in the preceding section in conjunction with a 
study on the conductivity of several distinct percolation models and has not only been 
applied to standard percolation. 

I compare the proposed in this paper TM algorithm(the modified algorithm) to the previ- 
ously published |lT6| , p!3| TM algorithms (the standard algorithms). As a standard algorithm 
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I used mostly the algorithm proposed in which was described in Sec II A. The code pub- 
lished in ||T6| was run with a technical improvement (zero elements check in the most-inner 
loop) only to be seen that the performance scaling is the same for both "standards". It is 
worth to note that the performance scaling ~ in 3D, reported in |]T6|, is an overestimation 
probably due to that technical item. 

In Fig. 1 1 display the amount of computer time required for the conductivity calculations 
by the modified and the standard algorithm on two different percolation models, namely or- 
dinary site percolation and one step bootstrap percolation In the bootstrap percolation 
model one generates a site configuration in several steps. First, one randomly occupies 



a specific small fraction of the lattice sites. Subsequently, one determines all empty lattice 
sites with at least two occupied neighbors and occupies these empty sites as well. The steps 
are repeated until no empty sites with two occupied neighbors remain. If such procedure 
stops after its first step I call it one step bootstrap percolation. The percolation-transport 
and correlation-length exponents of the one step bootstrap percolation model almost equal 



those of ordinary site percolation |23|. The computer time for these two models, when run- 
ning the standard algorithm in 3D, scaled in the same way - even with the same prefactor, 
so the averaged results are given on one curve ("3d") on Fig. 1. This coincidence encouraged 
using the data available in 2D for the comparison in the next paragraph. ( In 2D the 



two algorithms were applied to different models: the standard algorithm to the bootstrap 
model and the modified to ordinary percolation. ) 

As expected from the arguments in the preceding section, the modified algorithm displays 
the better scaling properties throughout. In 3D site percolation the computer time scales 
as ~ Lr"-"^^ for the standard algorithm, whereas the modification needs time proportional 
to L^'^° only. Similarly, in 2D we observe that the standard requires time ~ L^ ''^ whereas 
~ L}'^^ suffices for the modification. One can see that these values are appreciably smaller 
than the respective upper bonds given in Sec. 3. 2. 

The errors in the above values, given by the LSQ fitting procedure were smaller than the 
uncertainty coming from the eventual correction-to-scaling terms. A more careful analysis 
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is needed, but a reasonable value for the error should be within 0.1 - 0.2. All the tests were 
made at the percolation threshold for the respective model. The strips(bars) length was of 
the order 10^ in two dimensions and 10^ in 3D. Several strips were calculated for each width 
and model. If one have to compare with statistics accumulated by means of other method on 
square(cubic) samples, a strip 10^ x 100 may correspond to several thousend runs on (say) 
100^2 X 100^2 sample. (In Ref. H was found numerically that at pc the average length 
of the spanning clusters is around 2.0L and their number is = (3/8) (M/L) where M is the 
strip length.) 

The modified model was applied as well in studying properties of some "percolation- 
generated" fractals. After the extraction of spanning clusters at the percolation thresholds 
one adds new sites on the cluster perimeter in order some aerogell structures to be modelled 



In Fig. 2 I display timing results for such fractals (cf., p3[). The data sets (I) and (II) 
correspond to two ensembles with higher (I) and lower (II) fraction of additionally occupied 
perimeter sites of the spanning clusters. The set (III) has been taken for a model in which one 
considers second nearest neighbors in the spanning clusters as connected. As seen depending 
on the model the computer time scaling may vary significantly. 



IV. CONCLUSION 

The modified TM algorithm proposed in this work reduces significantly the computa- 
tional efforts required for obtaining the conductivity scaling for fractal structures in 2D 
and especially in 3D. In contrast to the L^ "^- computer time requirements of conventional 
TM algorithms in 3D the time requirements of the algorithm proposed in this work scales 
approximately as L^'^ for percolation clusters at pc- 

Extracting the percolation cluster backbone, instead of spanning clusters only, would 
ensure a further improvement of the performance. 

Probably the main disadvantage of the modification of the TM approach described here 
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is the complexity of the algorithm. Therefore, I have made the program publicly available 
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Figure captions 

Fig. 1 Comparison between the performance with and without extracting spanning 
clusters. (Time units are Sun SPARC 10 workstation CPU seconds) 

Fig. 2 Performance scaling for modified percolation models where, after extraction of 
spanning clusters, additional loops were closed: (1) by addition of new sites on the cluster 
perimeter; (II) by adding sites as in (I) but with a lower density, and (III) by increasing the 
connectivity range to second neighbors. (Time units are Sun SPARC 10 workstation CPU 
seconds.) 
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